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We investigate a gas of wet granular particles, covered by a thin liquid film. The dynamic evolution 
is governed by two-particle interactions, which are mainly due to interfacial forces in contrast to dry 
granular gases. When two wet grains collide, a capillary bridge is formed and stays intact up to a 
certain distance of withdrawal when the bridge ruptures, dissipating a fixed amount of energy. A 
freely cooling system is shown to undergo a nonequillibrium dynamic phase transition from a state 
with mainly single particles and fast cooling to a state with growing aggregates, such that bridge 
rupture becomes a rare event and cooling is slow. In the early stage of cluster growth, aggregation 
is a self-similar process with a fractal dimension of the aggregates approximately equal to Df ~ 2. 
At later times, a percolating cluster is observed which ultimately absorbs all the particles. The final 
cluster is compact on large length scales, but fractal with Df ~ 2 on small length scales. 

PACS numbers: 45.70.-n, 47.57.-s, 61.43.Hv 
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I. INTRODUCTION 

Granular materials are systems of macroscopic parti- 
cles which interact only when they are in mutual contact, 
and the interaction is dissipative. In spite of this simple 
definition, collective phenomena arising in such sytems 
are of utmost complexity, and have inspired strongly in- 
creasing research activities in recent years. The partic- 
ular interest in granular systems is mainly due to the 
fact that their importance spans from technology and ap- 
plied research to very fundamental questions of interdis- 
ciplinary relevance. On the one hand, storage and han- 
dling of bulk solids is among the most significant tasks in 
industrial technology, and still poses a large number of 
unsolved problems |l|, H, 0] ■ On the other hand, granular 
systems provide a comparatively simple, experimentally 
accessible model for physics far from equilibrium f3iH@|- 
This is at the heart of self-organization and pattern for- 
mation processes, so that granular systems have been 
considered as genuine model systems for structure for- 
mation on various length scales, including the formation 
of planetesimals from interstellar dust and the formation 
of planets and stars from accretion discs [7| . 

In most studies so far, models were inspired by dry 
granular systems, where the dissipative contact interac- 
tion consists in the loss of a certain fraction of the ki- 
netic energy in every impact. Adding a small amount 
of liquid to the granular system changes its properties 
dramatically: while dry sand can flow freely similar to a 
liquid, wet sand has properties of a plastic solid. This 
difference in the macroscopic behavior is reflected in a 
corresponding difference in particle interactions ^J- The 
collisions of dry granulates are typically purely repulsive 
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and characterized by the coefficient of restitution e which 
specifies which fraction of the kinetic energy is dissipated. 
Wet granular particles are covered by a thin liquid film. 
When two particles come into contact, the films merge 
and a capillary bridge is formed, exerting an attractive 
force on the particles. As the particles separate from 
each other again, the bridge stays intact up to a critical 
distance dc- At this point the bridge ruptures ^ and a 
fixed amount of energy is dissipated. Thus wet granu- 
lar particles are characterized by a hysteretic attractive 
interaction and a well defined energy which is dissipated 
when a capillary bridge ruptures. 

The existence of a well defined energy scale (and cor- 
responding time scale), which is absent in dry materials, 
is the essential microscopic ingredient not only of wet 
granulates but also of cohesive gases. In fact the liquid 
bridge can be thought of as a particular realisation of 
a more general cohesive force. A particularly important 
aspect of free cooling in cohesive gases is the aggrega- 
tion process which sets in, when the kinetic energy falls 
below the bond breaking energy. Wet granular systems 
may provide a realisation of various aggregation models 
and so-called sticky gases where particles move dif- 
fusively or ballistically until they collide and get stuck to 
an aggregate which is thereby growing. Such models have 
attracted a lot of interest [H El, M, M, Q d, d, [13, d, 
IT9I I , due to a wide range of applications ranging from the 
formation of dust filaments, snowflakes and clouds to the 
size distribution and impact probability of planetasimals 
in accretion discs. 

Kinetic properties of granular gases have been dis- 
cussed mainly for dry materials. In particular, free cool- 
ing has been studied extensively [l^, [2l|, and it was 
shown that the dissipative interactions are responsible 
for many novel phenomena, unexpected from the kinetic 
theory of molecular gases: The particles' velocities are 
not distributed according to a Maxwell- Boltzmann dis- 
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tribution '2'^, equipartition does not hold |23l. \24. \2l 
spatially homogeneous state is generically unstable [2f 
and linear and angular motion are correlated [2^ . 

Much less is known about wet granular media, which 
have been addressed only recently [8,, 2^, US H^j EH, IH, 

t ls^ . focussing on nonequilibrium phase transitions 
, the equation of state 33] , agglomeration ^2§, ^2^, , 
shear flow ^] , and cooling in one dimension [31, j32,] . 

Structure formation in wet granulates during free cool- 
ing has hardly been studied yet and is the focus of our 
paper which is organized as follows. In Sec. [iTl we in- 
troduce the model and discuss the decay of the average 
kinetic energy in Sec. IIIII Aggregation is discussed in 
Sec. IIVI before we present conclusions in Sec. |Vl A short 
summary of our results has appeared in [36| . 



II. MODELS 

In the present article, we are interested in the zero- 
gravity free cooling dynamics of wet granular gases. We 
assume the particles to be covered by a thin liquid film, 
as it is the case if the liquid completely wets the par- 
ticle material [s^]- The particles approach freely, until 
these surface films come into contact. The liquid then 
rapidly accumulates around the contact due to the in- 
terfacial forces. A capillary bridge forms at the contact, 
exerting an attractive force on the grains due to its nega- 
tive Laplace pressure. This liquid bridge is stretched but 
stays intact (or even continues to grow) as the particles 
move apart. The attractive force thus remains until a 
certain critical separation dc is reached, where the liq- 
uid neck becomes unstable and ruptures. As mentioned 
above, the hysteretic formation and rupture of the bridge 
gives rise to a characteristic loss of energy, AE, which de- 
pends upon the thickness of the liquid film wetting the 
grains. 

In order to design a suitable model, a few words on 
the details of this process are in order. The formation 
of capillary bridges is quite fast in real systems. Be- 
tween typical grains of one millimeter diameter it takes 
less than a millisecond. It is clear, however, that this 
formation cannot in general be considered instantaneous 
if the velocity of the impacting grains, Vi, is large. If the 
time scale of the impact process, which may be written 
as dc/vi, is of the same order or even smaller than the 
time of capillary bridge formation, the accumulated liq- 
uid volume of the bridge, and hence AE, will be smaller 
than for slow impacts. However, this will not greatly af- 
fect the main features of the wet system, in particular 
as to its characteristic difference from the dry granulate. 
In order to see that, we compare the effective restitu- 
tion coefficient of the dry and of the wet system. This 
is shown in Fig. [Tl where the restitution coefficient for 
the dry system is shown as the dotted curve. It tends 
to be mildly depending on impact energy [4], Ei, with a 
negative slope throughout. The e ffective restit ution co- 
efficient of the wet system, Ceff = — AE/Ei, is shown 



as the solid curve, assuming constant AE. In strong con- 
trast to the dry system, it has a zero at AE/E[ = 1, and 
a markedly positive slope. This illustrates the dramatic 
difference between these two systems. The dashed line 
qualitatively accounts for the effect of finite formation 
time of the capillary bridge. Since £ofr must stay be- 
low one, the difference between the solid and the dashed 
curve is very limited, and the qualitative picture concern- 
ing the comparison of dry and wet granular gases remains 
unchanged. 
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FIG. 1: Restitution coefficients for dry (dotted) and wet (solid 
and dashed) granular systems, plotted vs the impact energy 
in units of the wet energy loss, AE. The main feature in 
the wet case is the zero at Ei = AE, which is unchanged if 
the finite formation time for capillary bridges is taken into 
account (dashed curve). 

Our system consists of N identical and spherical parti- 
cles with diameter d and mass m in a three-dimensional 
cubic volume V = L^. The particles have a hard core 
interaction, such that two particles are reflected elasti- 
cally, if their centers of mass reach the hard-core distance, 
which is the particle diameter d. 

To account for the liquid film, a liquid bridge is allowed 
to form between a pair of particles if they come close 
enough ("close enough" is specified later). When these 
particles are moving apart and their distance exceeds the 
bond breaking distance dc, the liquid bridge will break 
and a fixed amount of kinetic energy AE is dissipated; 
thereby, momentum is conserved and the relative velocity 
Vrci changes to w^^j according to 



AE 



(1) 



with the reduced mass /i = m/2. If, however, the rela- 
tive kinetic energy is smaller than AE', the particles are 
elastically reflected towards each other. The effect of the 
capillary force, which is present in reality for distances 
up to dc, is thus solely modelled by the enrgy loss which 
occurs when d = dc. This has been shown before to be a 
very good approximation (33 |. and enables event-driven 
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simulations as discussed below. For the formation of the 
liquid bridge, we distinguish between two models: 

In the thin film model, the liquid bridge forms when the 
particles touch, i.e. the distance of their centers is equal 
to d. This model assumes that the liquid film covering the 
particles is infinitesimally thin and the capillary bridges 
form a thin liquid neck, which breaks off at the critical 
distance d^. 

In the thick film model, a liquid bridge forms as soon 
as particles come closer than the critical bond breaking 
distance c?c- This model assumes that the outer diameter 
of the liquid film is dc and its shape stays spherical and is 
not deformed by the particles. Although this may seem 
unphysical, we include this case in our study because 
similar assumptions have been used in many simulation 
studies in earlier articles. As it will turn out, the differ- 
ences in most of the results are only minute. The two 
models are illustrated in Fig. [51 
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simulations of the wet system with full molecular dynam- 
ics simulations integrating the equation of motion [s^ . 
We found good quantitative agreement in the results of 
both methods, justifying the event-driven approach we 
chose exclusively for the present study. 

We use dimensionless units such that AE = 1, par- 
ticle mass m = 1 and particle diameter d = A. The 
bond-breaking distance is chosen as dc — l.OTd, unless 
noted otherwise, and volume fraction, — ird"^ /Q ■ N/V, 
is varied from « 0.06% up to 15.6%. We use periodic 
boundary conditions in the x- and y-direction and hard 
walls in z-direction. 



III. COOLING DYNAMICS 

We define the granular temperature T — jj^ 12iLi fnwf 
and investigate its decay in time from a given initial value 
To ^ AE. In all our simulations we choose Tq = 45Ai?. 
Simple arguments can be used to derive an analytical 
form of the temperature decay. In each coUsion a cap- 
illary bridge ruptures with probability Pbb, giving rise 
to dissipation of a fixed amount of energy, the bond- 
breaking energy AE. Particles collide with frequency 
/coll, so that the average loss of energy per unit time is 
given by: 

-— -/coll- A£;-n&. (2) 

2 dt 2 ^ ' 

The factor i takes into account that two particles are 
involved in one bond-rupture. 



FIG. 2: Illustration of the thin film model and thick film 
model. In the thick film model, the liquid bridge forms, as 
soon as the bond breaking distances dc overlap. The same 
initial configuration in the thin film model does not create 
a liquid bridge, since the hard cores of the particles do not 
touch. Thus, the particles just pass by. 

In general there is some energy being transfered to the 
atomic degrees of freedom of wet grains as well. In this 
paper we are going to neglect this dissipation mechanism 
because it is usually small as compared to the energy 
loss due to the breaking of capillary bridges, especially if 
the granular temperature is small. However, we want to 
point out that such a dissipation mechanism can easily 
be incorporated in the simulations, replacing the elastic 
reflection by incomplete normal restitution. We restrict 
ourselves here to perfectly smooth particles, such that 
translational and rotational motion are decoupled. Fur- 
thermore, we investigate free cooling only, so no gravity 
is present, and no energy is injected into the system. 

The particular way of accounting for the liquid film 
used in these models makes it possible to use an event- 
driven simulation scheme. The possible events are the 
reflection of the particles at the hard core distance d and 
the crossing of the bond-breaking distance d^. As men- 
tioned above, we have previously compared event-driven 



A. Early stage of cooling 

In the early stage of cooling the average kinetic en- 
erghy per particle is much larger than the bond breaking 
energy, so that Phf, « 1 and almost every collision gives 
rise to dissipation by AE. For a dilute gas, the colhsion 
frequency 



/coll = Agid)an 



T 



(3) 



is well established, with the particle density n — N/V 
and the pair correlation function at contact g(d) — 
J"(lZ'!p)3 (6-g- [4])- The two models differ only in the cross 
section a (see Fig. [2]), which is given by a = d'^n in the 
thin film model and a = d^ir in the thick film model. 

The only temperature dependent quantity remaining 
on the right hand side of Eq. ^ is the collision frequency, 
/coll ^ VT from ([3]), giving rise to the following simple 
equation: 



dT 
It 



-VT, 



(4) 



which is solved by T{t) ~ (t— tg)^- Insertig the prefactors 
and the initial value T(0) = Tq, one obtains, similar to 
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Haff's law [38j, an analytical form of the decay of the 
temperature: 



To (1 - t/to)^ for t < to 







for t > to 



with a charecteristic time scale 



to 



2g{d)anAE 



(5) 



(6) 



Note that, in this simplified model, the assumption 
that every collision causes an energy loss AE gives rise to 
a time-scale to after which all energy is dissipated. Even 
though this assumption does not hold for all times in the 
simulation (since the bonds do not break anymore if the 
relative kinetic energy is too small), the timesacle to has 
a clear physical relevance. It sets the time after which 
the temperature is comparable to the bond-breaking en- 
ergy AE and after which persistent clusters will form. In 
Fig. [3] the evolution of the granular temperature T from 
the simulation is compared to ([5]) for different volume 
fractions: 0.061% < ch < 15.6%. 




FIG. 3: (color online) Decay of the granular temperature T for 
the thick film model and volume fractions (from left to right) 
(/) = 15.6%, 7.81%, 3.90%, 1.95%, 0.98%, 0.49%, 0.24%, 0.12%, 
0.061%. N — 262144 particles are fixed. The corresponding 
solid lines show the analytic form (O with a decay to zero at 
time to, given in (|6l). At that time the temperature of the 
simulated granulate shows a rapid transition to a value below 
the bond- breaking energy AE = 1. In the inset temperature 
data are plotted versus scaled time t/to, such that data for 
different volume fractions collapse onto a single curve. 

In the simplified cooling law ([5]), the volume fraction 
only enters into to- Hence we try to superimpose the data 
by scaling time with to- As can be seen in the inset of 
Fig. [3l the data obey the expected scaling well, except 
for the long time limit, which has different asymptotic 
behavior and is treated in the next section. 

In Fig. m we compare data from the thin and thick film 
model for two volume fractions, (j) = 1.95% and 0.24%. 
The difference is solely due to different scattering cross- 
sections, entering in to ^ and can be absorbed into the 
rescaling of time by to- 
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FIG. 4: (color online) Decay of the granular temperature T 
for the thick film model (•) and the thin film model (o) for 
volume fractions 6 = 1.95% and 6 — 0.24%. 



B. Late stage of cooling 

In the late stage of aggregation, when the system is 
strongly aggregated, it becomes very unlikely that a cap- 
illary bridge ruptures. Hence we observe a very slow time 
evolution of our system. The slow decrease of the tem- 
perature can be understood with simple arguments. The 
probability Pbt to break a bond is given by the probabil- 
ity to find a kinetic energy larger than AE: 



Pi 



bb 



d^v eimv^ /2 ~ AE) w{v). 



(7) 



We approximate the velocity distribution w{v) by a 
Maxwellian 



w(v) 



\2TTT{t) 



3/2 



,-,nv^/(2T(t)) 



(8) 



and evaluate the above integral in the limit T{t)/AE 
0. The probability to break a bond becomes exponen- 
tially small in that limit: 



P, 



bb 



1/2 



^-AE/T 



(9) 



The decrease of kinetic energy, as given by Eq. ((2), is 
now dominated by the probability to break a bond. The 
collision frequency /cou is not known for the clustered 
state, but is expected to be proportional to T^/^. Using 
© and /coU oc T^/^ in the rate equation ^ yields: 



dT 
'dt 



- 76 



-AE/T 



(10) 



The prefactor 7 is determined by the precise form of the 
collision frequency. Separation of variables can be used 
to integrate Eg. pO]) 



(11) 



j-T/AE 

/ dxe^/^ = -7(i-ti) 

Jt, /AE 



iTi/AE 
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with the intital value Ti = T{ti). In the asymptotic 
Hmit r ^ and Ti with T < Ti, one finds a 
logarithmically slow time decay of the temperature 



T 

AE 



1 



Ht) 



(12) 



which is due to the very low probability to break a bond, 
Eq. This is in strong contrast to the algebraic time 
decay observed for dry granular systems with coefficient 
of restitution e < 1. [J 

In Fig. [SI the full solution (fTT|) is compared to the 
simulation data, showing good agreement. The unknown 
prefactor 7 is a fit parameter. It is noteworthy that for all 
densities, the temperature seems to approach a universal 
curve as t 00. 
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FIG. 5: (color online) Asymptotic time dependence for several 
volume fractions as in Fig.O data (dots) in comparison to the 
analytical results (lines) 



C. Partitioning of the energy into translational, 
rotational and internal degrees of freedom 

After the time to has passed, stable clusters emerge. 
For the definition of a cluster, we define particles as 
neighbors, if a bridge is formed and the relative kinetic 
energy is not sufficient to break it. This makes sure that 
particles which are just "passing by" , are not consid- 
ered neighbors. A cluster is a set of particles connected 
through this neighbor-relationship. Hereby we refer to 
the cluster mass m as the number of particles a cluster 
contains. Clusters defined in this way are not truly sta- 
ble. Particles belonging to the cluster are occasionally 
kicked out, if hit by a very energetic particle. 

For a more detailed understanding of the system, we 
investigate the cooling dynamics on the cluster level, and 
determine how energy is partitioned among the degrees 
of freedom. We split the total temperature T into three 
constituents, the translational temperature defined via 
the center-of-mass velocities of the clusters, the rota- 
tional termperature defined via the angular momenta of 



the clusters, and the internal temperature describing the 
relative movement of the particles inside a cluster. These 
three temperatures are defined as follows. 

Our definition of neighborhood relations gives rise to 
Uci distinct clusters numbered by i = 1, rici- We denote 
by Mi the i-th cluster with particles. Its centre of mass 
position and velocity are given by: 



rui ^ — ' 



and V, 



(13) 



Note that single particles with nii — 1 are also considered 
as clusters. 

The center of mass movement of each cluster has 
/trans, i — 3 translatioual degrees of freedom, so that the 
total number of translational degrees of freedom of these 
clusters is simply 3nci. Homogeneous cluster translations 
are thus characterized by the translational temperature 



Ttra 



2 

3nci 



(14) 



Analogously, the rotational temperature describes the 
energy in homogeneous cluster rotations. The angular 
momentum, L.^, of cluster i is given in terms of the rel- 
ative particle positions f = R.^ — and velocities 



(15) 



The rotational energy of cluster Afi with > 2 is thus 
given by 



Erot.i — nLil. "'^ Li 



(16) 



where the moment of inertia tensor I. is defined in the 

—I 

usual way. The case mi = 2, requires special treatment, 
since the inertia tensor is singular. The rotational energy 
of a dimer can be easily calculted to i?rot,i = (vi— v2)5^/4, 
where (vi —V2)± denotes the relative velocity perpendic- 
ular to the axis of the dimer. The rotational temperature 
is thus 



Trot ■ — 



(17) 



with /rot.i — 2 for dimcrs and /rot,i — 3 for larger clusters. 

All the left-over kinetic energy Eint describes the rel- 
ative movement of particles inside a cluster and con- 
tributes to the internal temperature. Each cluster has 
a total of 3mi degrees of freedom, so that the remain- 
ing number for internal degrees of freedom is /int,i = 
3TOi - /trans.i " /rot.i- The internal temperature Tint is: 



Ti,- 



int. 2 



(18) 
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FIG. 6: (color online) Top: Division of the total 3N degrees of 
freedom into the translational, rotational and internal parts, 
dependent on time. Bottom: Evolution of the total (•, black), 
translational (A, green), rotational (■, red), and internal (♦, 
blue) granular temperatures. Data for A*' — 262144 particles 
and volume fraction (p = 1.95% are shown; the behavior is 
qualitatively the same for all investigated system sizes. The 
horizontal line at 2/3 corresponds to the bond breaking en- 
ergy. 



Fig. [6] (top) shows how the total of 3N degrees of free- 
dom divide up into translational, rotational, and internal 
degrees of freedom. The corresponding temperatm'es are 
shown in the lower half of the figure. As one might ex- 
pect, for t ^ to almost all degrees of freedom are trans- 
lational, since most clusters are just single particles, and 
Ttrans ~ T. Keeping in mind that two particles are only 
defined as neighbors if their relative velocity is not suf- 
ficient to break the bond, only stable clusters (mostly 
dimers) enter the internal and rotational temperatures, 
and therefore T,ot,Tir,t < fAi; = | for </<o < 1. (ioj 

In the transitional regime t ~ to, when the num- 
ber of intermediate size clusters increases, the rotational 
degrees of freedom become important. Larger objects 
can have higher rotational energies without rupture [50| , 
therefore the growing clusters obtain energy from caught 
particles, and thus Tiot increases until reaching the value 
of Ttrans- After that, the energy of the incoming lumps 
is not sufficient to increase Trot any further. 

In contrast to the homogeneous cluster rotations, the 




FIG. 7: (color online) Snapshot of the system with volume 
fraction — 0.48% and A'^ = 262144 particles taken at time 
t ~ 12to; the largest cluster (grey) contains 22% of the parti- 
cles. Particles of the same cluster have the same color shade. 



internal degrees of freedom which have higher energies 
than AE will in most cases result in a bond rupture, 
independent of the cluster size. Therefore, Tint decreases 
monotonically. At late times t ^ to, large clusters have 
formed, thus almost all degrees of freedom are internal 
and T « Tint- 



IV. AGGREGATION 

When the average kinetic energy per particle is com- 
parable to the bond breaking energy, t tg, the sys- 
tem starts to form aggregates, which seem to grow in 
a self-similar process. In the following we are going to 
analyze these aggregates and compare them to cluster- 
cluster aggregation [33| models. As time proceeds, larger 
and larger clusters are formed. We observe a spanning or 
percolating cluster for all finite densities, and ultimately 
all particles and clusters have merged into a single clus- 
ter. 

Figs. [7] and [5] show snapshots of a system a.t t ~ 12to 
and t — 52to with small volume fraction, ip = 0.48%. 
At the smaller time the system is not yet percolating, 
even though rather large clusters have already formed, 
the largest one (in grey) contains 22% of all particles. 
The second snapshot, taken at a much longer time, shows 
a spanning cluster. At such large times the average ki- 
netic energy is much smaller than the bond breaking en- 
ergy (T « 0.06Ai?), so that bonds almost never break 
up. The cluster shown is already well beyond the criti- 
cal time for percolation with 99% of the particles in the 
cluster. 

Fig. [5] shows the evolution of the cluster mass distri- 
bution Nm (t) , which is the number of clusters containing 
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FIG. 8: (color online) Same as Fig.[7]for t ~ 52to; the largest 
cluster contains 99% of the particles. 
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FIG. 9: Histogram of the cluster mass distribution dependent 
on time, for volume fraction (j> = 3.9% and = 262144. The 
number of clusters at the respective time and size is color 
coded on a logarithmic scale so that the single largest cluster 
is visible. At t « 2.5to one can see the large cluster emerging, 
clearly distinguishable from the rest of the distribution. 



m particles at time t. One can clearly see that after 
some time, t « 2Mq, which depends on volume fraction, 
the largest cluster emerges from the rest of the distribu- 
tion. For all volume fractions a gelation transition was 
observed at the percolation time > to- The critical 
behavior of the gelation transition is still controversial. 
Since aggregation is a nonequilibrium process, there is a 
priori no reason that it should be in the same univer- 
sality class as the corresponding equilibrium percolation 
transition. Yet there is some evidence in favour of this 
conjecture. Gimel et al. [40| observe a crossover from 
self-similar growth at small times and volume fractions - 
called the flocculation regime - to the percolation regime 
around tc- In the latter they observe critical exponents as 
in standard percolation theory. Kolb and Herrmann 41] 
on the other hand obtain values for the fractal dimen- 
sion of the percolating cluster, distinct from percolation 
theory as well as from flocculation theory. Both studies 
refer to diffusion limited cluster-cluster aggregation. 

In this paper we do not analyze the gelation transi- 
tion in detail but defer such a discussion to future work. 
Instead we investigate two regimes in detail in the fol- 
lowing: 

a) The self-similar growth process, or flocculation 
regime, which is present for small times and volume frac- 
tions. 

b) The properties of the final cluster which emerges, 
when (almost) all particles have aggregated to form one 
large cluster. 



A. Self-similar growth 

1. Fractal dimension of the aggregates 

A central quantity of aggregation models is the frac- 
tal dimension of the aggregates. It is usually deter- 
mined from the radius of gyration as a function of cluster 
mass. We consider a cluster of m particles with posi- 
tions (ri, Tm) and define its radius of gyration by (see 
e.g. 0) 

^ m ^ m 

rg(m) = — y^(ri - f)^ with f = — . (19) 

i=l i=l 

If the clusters are fractal we expect a scaling relation for 
large m of the form 

rg - m^/^f (20) 

which yields the fractal dimension Df. This method is 
commonly used in aggregation models, where particles 
move diffusively, ballistically, or are interact ing and stick 
to the aggregate once they touch it [13, [H, [l3, SI] . 

In Fig. [To] we show the radius of gyration for a sytem of 
262144 particles at volume fraction (j) = 1.96%. Several 
snapshots of the ensemble of growing clusters have been 
taken at times to < t < tc with the percolation time tc, 
when a spanning cluster is first observed. The data scale 
well according to Eq. ((20|l . some scatter is observed for 
the largest masses, corresponding to times close to the 
percolation transition. 

In contrast to aggregation models, where the clusters 
are static and do not break up, we occasionally do ob- 
serve the breaking of bonds. In addition there are inter- 
nal deformations of the clusters during growth, so that 
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FIG. 10: (color online) Radius ol gyration as a lunction ol 
cluster size for a system of 262144 particles at volume fraction 
(/> = 1.96%; different colors/shades correspond to simualtion 
times between to (yellow/light gray) and 4io < tc (black)); 
The slope of the solid line corresponds to Df = 2; inset: fractal 
dimension as a function of time, extracted from the slope of 
the curves in the main figure. 



the fractal dimension could depend on time. We have 
therefore checked the relation between m and Tg{m) for 
many instances of time and show the fractal dimension 
as a function of time in the inset of Fig. [TOl As can be 
seen from the Figure, there is no systematic dependence 
on time, and the fractal dimension is close to D f = 2. 



2. Cluster size distribution 

All information about the connectivity of the clusters 
is contained in the cluster size distribution Nm{t), the 
number of clusters of size m at time t. In Fig.[Tl]we show 
Njn{t) for a system with (j) = 1.96% and N = 1048576. 
The time interval has been chosen such that to < t < 
2to < tc ~ 4^0 (for this volume fraction). In this time 
interval the mean cluster mass increases roughly by a 
factor of 30. 

It has been suggested [e.g. 43]) that for aggregat- 
ing systems the mass distribution evolves towards a self- 
preserving scaling form, independent of the initial distri- 
bution: 



Nm{t)=m-^f{m/m{t)) 



(21) 



where the time dependence is only contained in the mean 
cluster mass 



(22) 



This scaling form has been applied sucessfuUy to various 
aggregating systems [l^, [H, [lE [3 IH, El, SE] , involving 
fractal as well as non-fractal objects. Mass conservation 
requires 6 — 2 [43]. 

We plot in Fig. [T^] the scaling function f{m/fh) = 
Nm{t)m'^ for the same data sets as in Fig. [11] We expect 



FIG. 11: (color online) The cluster mass distribution N,n{t)- 
The different graphs represent different times, which are in- 
creasing from top to bottom (left side of the graph). The 
inset shows how the mean cluster mass increases during the 
investigated time period. 



scaling to hold only in the aggregation regime, i.e. for 
times not too close to tc, where the system gels (see 
sec. IIV A3I) . Hence we restrict ourselves in Fig. [T^ to 
times to < t < 2to. We have also left out the data 
points for m = 1, i.e. clusters consisting of single par- 
ticles. As can be seen from Fig. [12] the data scale very 
well. Deviations occur only for times close to the perco- 
lation transition (not shown here), where they should be 
expected. 
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FIG. 12: (color online) Rescaled cluster size distribution 
f{m/rh) = Nm{t) ■ from eq. (|21|l versus the normalized 
cluster mass m/fh. The color coding as in Fig. Illl is used. 



3. Number of clusters 

Another characteristic of a realisation of clusters 
is simply the total number of clusters nci(t) — 
X]m=i-^™(^)' which decreases as aggregation proceeds. 
As long as the system is in the scaling regime (i.e. rela- 
tion ([^T]l is fullfilled), the mean cluster mass, in{t) and 
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the number of clusters are simply related: rh(t) ~ n~i^. 
However, as mentioned above, the scaling relation (pTj) 
only holds in the aggregation regime and is expected to 
break down as the percolation transition is approached. 
At that point, to should diverge due to the formation 
of a spanning cluster. On the other hand, there is still 
a large number of smaller clusters coexisting with the 
macroscopic cluster, so that Ud/N remains finite at the 
percolation transition. 

The aggregation of particles to larger objects has 
been investigated for various ballistic aggregation models 
[T3I [Tsl . [46j , where spherical particles of mass m — 1 
and diameter d = do move ballistically, until two of them 
collide to form clusters irreversibly. In a particularly 
simple model, one assumes that two colliding particles 
form one larger spherical particle with conserved momen- 
tum and a mass to equal to the sum of the two particles 
masses, so that m is always equal to the number of initial 
particles contained in a given cluster. For spatial dimen- 
sion D, the diameter increases like d = m}/^do, assuming 
the particles to be compact spheres which conserve vol- 
ume when merging. For this model^a mean field theory 
and simple scaling arguments [l^, [ll| yield the de- 
pendence of the expected average mass to on time like 
TO ^ with an exponent ^ = 2D/{D + 2) (assuming 
to = 0). 

Since the aggregating clusters in our system are not 
compact, but fractal objects with fractal dimension Df, 
the assumption for the diameter d ~ rri^l^ does not hold 
and must be changed to c? ~ to^/^' . With this assump- 
tion, we follow the scaling arguments of Trizac ct al. [l^| , 
and find the scaling relation between to and t. 

We assume that the number of clusters per volume, 
rici, is reduced by one whenever two clusters collide: 



dnc\/dt -/coll • nci 



(23) 



The collision frequency [4| is approximately given by 
/coll ~ d^~^nc\v with d oc rg the linear dimension of 



the cluster and v its typical velocity. The average mo- 
mentum should scale as p ^ m^l"^ |18| . and therefore 



V — p/j 



^1/2 



1/2 



Plugging in all these scaling relations as well as to • 
one obtains: 



"dT 



c\- V - u r 

which is solved by 

ncl~ (t-t*)"''''/<''''-"'''+'\ 



5/2-(_D-l)/Df 



(24) 

r-Dt 
g ' 

(25) 
(26) 



where the integration constant t* is the onset of cluster 
growth. In our context t* « to This implies the 

following growth law for the mean cluster mass in the 
scaling regime: 



TO 



(t-t*)^ with e 



2Df 



3Di-2D + 2' 



(27) 



which generalises the result for compact objects, ^ = 
2D/{D + 2) with D = D{ to fractal ones with D ^ Df. 

In Fig. [13] we show how the number of clusters de- 
creases over time as larger and larger aggregates form for 
t > to. 




FIG. 13: (color online) Evolution of the mean cluster mass. 
Labeling and parameters as in Fig.|3l For the inset, the origin 
of the time-axis has been shifted to the transition point t* to 
investigate the scaling relation rid ~ {t — t*)~^. The solid line 
has a slope of —2. 

The inset of fig. [13] investigates the scaling behavior 
(PS)) . with the origin of the time axis shifted to the tran- 
sition point t*. One can see that the slope of ^ = 2, 
obtained from ((27)) for D = 3 and _Df = 2 is in good 
agreement with the simulation. 



B. Properties of the asymptotic cluster 

The fractal dimension of the largest cluster - well be- 
yond the percolation transition for most volume fractions 
- will be the main focus of this section. In particular we 
determine its fractal dimensions and coordination num- 
bers. 



1. Fractal dimension from radius of gyration 

One way to determine the fractal dimension is the ra- 
dius of gyration, as was done in Sec. lIV ATl for aggregates. 
Here, however, we only have one large cluster and have 
to find a way to obtain the function rg (to) as a function 
of cluster size to. We implement this in following way: 
Starting from a random particle of the cluster, we mark 
all particles that can be reached through i neighbor-to- 
neighbor steps. Thus, for every i, we get a partial cluster 
with m{i) particles and radius of gyration rg(i), which 
yields the scaling relation rg ^ to^/-^*' and the fractal 
dimension D{. For good statistics, we repeat this proce- 
dure 100 times (each with a different initial particle) and 
average over the obtained values of Vg. Note furthermore 
that the procedure takes care that no particle is marked 
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FIG. 14: (color online) Radius of gyration dependent 
on the mass of the partial cluster at simulation time 
t 27to- The particle number is fixed = 262144 
and the volume fractions are (from bottom to top) = 
15.6%, 7.81%, 3.90%, 1.95%, 0.98%, 0.49%. The lines along 
the data points are the respective fits. The outer solid lines 
have slopes 1 /2 (top) and 1 /3 (bottom) corresponding to frac- 
tal dimensions of 2 and 3, respectively. 



a second time, in order to make sure that one does not go 
through the cluster several times because of the periodic 
boundary conditions. 

In Fig. [13] we show the results of this procedure for 
the radius of gyration rg as a function of m for different 
densities. For high volume fractions we are well beyond 
the percolation transition and hence expect D[ — 3 on 
the largest length scales of the cluster. This is clearly 
seen in Fig. [Til e.g. for 4> — 15.6% and 10^ < m < 
10^. On smaller length scales, however, we find a fractal 
dimension Df « 2. For smaller volume fractions, the 
crossover to Df = 3 happens at larger masses and hence 
the "interior" region extends to larger scales. 



2. Fractal dimension from box counting algorithm 

To further investigate the Hausdorff dimension of the 
largest cluster at intermediate length scales, we use the 
box counting algorithm [43, 13 ■ The system is divided 
into sub-boxes of edge length Lbox- Then each box which 
contains or hits at least one particle is marked. In this 
way, we find the number of boxes iVbox necessary to cover 
the whole cluster. This number should scale with Lbox 
like 

iVbox ^ ibox'' - (28) 

with the Hausdorff dimension D{. 

On length scales much smaller than the particle di- 
ameter, Lbox ^ d, the system obviously behaves three- 
dimensionally. In this regime, the number of filled boxes 
-^box is just the volume fraction times the total number 



of boxes iVbox,tot = i^/^box' therefore: 

iVbox = . (29) 

Since our system is finite and contains a system- 
spanning cluster, the scaling behavior on large length 
scales Lbox ~ L should also be three dimensional. On 
this length scale, almost all the boxes should be filled, so 
that 

A^box = 7^ . (30) 

-'^box 

In particlar, the relation must include the point 
(Lboxj^box) = (L, 1), since a box of the system size in- 
cludes all particles and will certainly be marked. 

Only in the regime between these two limiting cases is 
it possible to observe the fractal dimension with the box- 
counting method. Comparing pS)) and shows that 
the interesting range is proportional to | log (f>\, which only 
depends on the volume fraction, but not on the particu- 
lar choice of the system size. A schematic plot is given 
in Fig. [Tni where the number A^box of boxes containing 
particles is plotted against the edge length Lbox of a box. 




FIG. 15: Schematic double logarithmic plot of the box size 
ibox versus the number of boxes A'^box of that size needed to 
cover the cluster. The negative slope is the fractal dimension. 
We expect three scaling regions: For small and large Lbox, the 
system should behave three dimensionally, and the region in 
between yields the non-trivial fractal dimension. If only the 
particle centers are considered, the algorithm simply counts 
the number of particles in the cluster for Lbox ^ d resulting 
in a horizontal line (dotted line). 

For numerical reasons, it is very tedious to observe the 
expected slope of —3 for small Lbox, because of the vast 
amount of boxes to account for. Since this regime is not 
relevant anyway, it has only been investigated exemplar- 
ily and is reached for Lbox ^ 0.03d. For all other runs 
we simplify the algorithm and only use the centers of the 
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particles, i.e. a box is only marked, if a particle center is 
inside. With this definition the number of boxes needed 
to cover the system for small box sizes Lbox < d is just 
the particle number N , resulting in a horizontal line on 
the left side of the graph, instead of the slope —3 (dotted 
line in Fig. [T5|) . 

1 5 10 50 100 



10"^ 




1 5 10 50 100 



FIG. 16: (color online) Top: A'^box versus Lbox at time t « 27to 
and volume fraction (j> = 1.95% for the box counting algo- 
rithm; particle number is varied from bottom to top, accord- 
ing to TV = 32768,65536,131072,262144,524288,1048576. 
The straight lines are fits to the data to the left and right 
of the cross-over point Leo ~ 25d, which is also a fitting pa- 
rameter and shown as a star (*). Bottom: As top, but A^box 
normalized by cluster mass; the solid lines have slopes —2 
(left) and —3 (right); the vertical dashed line represents the 
particle size. 

Fig. [16] (top) shows the outcome of the box-counting 
algorithm, at a time t « 27to, where roughly all parti- 
cles are inside the largest cluster. It yields the relation 
between the box size Lbox and the number of boxes of 
that size, needed to cover the cluster. The slope of that 
curve is the negative fractal dimension. The result for 
different system sizes, but with the same volume fraction 
are presented. As proposed, for all system sizes, there is 
a cross-over point Leo, at which the slope changes. On 
length scales between d and Leo, the fractal dimension is 
roughly 2 (the fits yield values between 1.92 and 2.03). 
Above Leo the fractal dimension has a trivial value of 
about 3 (fit values between 2.95 and 3.00), which means 
that on these large length scales all the boxes are filled 
and is therefore an indication that the cluster is system- 



spanning. 

In the lower half of Fig. (TH] we show the number of 
boxes normalized by the cluster mass. The data collapse 
well onto a single curve, obviously with the same slopes. 
Here one can see very well that for systems with the same 
volume fraction, the slopes as well as the cross-over point 
do not depend on the absolute system size. 

Results of the box counting algorithm for different den- 
sities are presented in Fig. [T71 We only include densities 
for which a spanning cluster has developed. The .slopes 
of —2 and —3 in the two scaling regions are not affected 
by the volume fraction 0, but the size of the non-trivial 
region (with D{ « 2) is seen to increase significantly as 
the density decreases. Even for the lowest density, the 
size of the scaling region is less than 2 decades, which 
makes it difficult to extract precise values for the fractal 
dimension. For the three most dense systems, the scaling 
region is less than one decade. As discussed earlier, this is 
an intrinsic feature of the "high" density systems, which 
can not be resolved by taking larger systems {N, L — > oo 
with constant (j))- As the inset shows, we can collapse all 
data on a single curve by rescaling iVbox with and 
Lbox with (j) in agreement with the dependence of the 
crossover length on | log0|. 
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FIG. 17: (color online) Result of the box counting algo- 
rithm, as in Fig. [in] at t « 45to; the particle number is fixed 
— 262144 and the volume fraction varies from left to right 
according to </> = 15.6%, 7.81%, 3.90%, 1.95%, 0.98%, 0.49%; 
the lines along the data points are the respective fits; the 
outer lines have slope —2 (top) and —3 (bottom). 



3. Coordination number 

Given the definition of a neighborhood relation (two 
particles are neighbors if they have built a bridge and 
their kinetic energy is not sufficient to break it), we can 
extract the average number of neighbors of a particle, 
i.e. the average coordination number. In Fig. [18] we show 
histograms for the coordination number in the percolat- 
ing cluster for two different bond breaking distances. As 
one would expect, these distributions are rather broad 
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with coordination numbers between one and thirteen. 
The smaller bond-breaking distance (left) gives rise to 
a more asymmetric distribution with more weight for 
smaller coordination numbers. 





FIG. 20: Influence of the bond-breaking distance dc on the 
final value of the average coordination number for the thin 
film model (right, (j) = 0.24% and N = 10648) and on the 
fractal dimension of the thick film model (left, system and 
colors as in Fig. I19|l . 



FIG. 18: Histogram of the coordination number for two differ- 

ent bond breaking distances dc = l.Old (left) and dc = 1.07d average coordination number^ut not influencing the 

(right); for both plots N = 262144 and = 1.95%. stucture on larger length scales [52]. 



In Fig. I19( we show the time evolution of the aver- 
age coordination number for different bond-breaking dis- 
tances dc- After a strong increase at the time to, the 
coordination number continues to grow slowly. This slow 
increase is strongly suppressed in the thin film model 
(right) as compared to the thick film model (left). Within 
the thin film model the slow growth with time is further 
suppressed for decreasing bond-breaking distance dc- 




FIG. 19: (color online) Time evolution of the average coordi- 
nation number of particles in the largest cluster (A^ = 262144 
and (j} — 1.95%); the critical break-off distances are dc — 
1.07d, 1.035d and l.Old from top to bottom; thick film model 
(left) in comparison to thin film model (right). 

As one can see in Figs.[THl[IH the coordination number 
becomes smaller for smaller dc. This is reasonable, be- 
cause the particles can more easily collect neighbors for 
higher dc- As dc ^ d the average coordination number of 
the thin film model approaches 6, which is the isostatic 
value. This is demonstated in the right half of Fig. [20l 
where we plot the asymptotic coordination number as a 
function of dc. Here the asymptotic value is taken, when 
T < OMAE for the first time. 

Naively one might expect that the increase of the co- 
ordination number with larger dc is caused by a com- 
pactification and therefore accompanied by an increase 
of the fractal dimension. However, as can be seen in the 
left part of Fig.[20l there is no significant influenece of dc 
on the development of the fractal dimension. Thus, we 
conclude that this compactification is mostly occuring on 
the single particle length scale and therefore increasing 



V. CONCLUSIONS 

We have analysed a simple model of a wet granulate 
allowing for large scale event driven simulations. A cen- 
tral feature of wet granulates is the existence of an energy 
scale AE associated with the rupture of a capillary bridge 
between two grains. This energy scale has important con- 
sequences not only for the phase diagram [sj but also for 
the free cooling dynamics investigated in this paper. The 
most important feature is a rather well defined transition 
at a time Iq, when the kinetic energy T of the particles 
becomes equal to AE. 

For t < to the particles are energetic enough to supply 
the bond breaking energy AE, so that very few collisions 
result in bound pairs and most particles are unbound. 
Cooling is very effective in this regime, but drastically 
different from a dry granulate. Whereas in dry granulates 
the dissipated energy is proportional to the energy of 
the colliding particles, in wet granulates the dissipated 
energy is AE, independent of the energy of the colliding 
particles so that T ^ \/T . Consequently Haff's law does 
not hold and is replaced by T{t) = T(0)(1 - tjtof for 
t < Iq. The simulations are in very good agreement with 
this cooling law for t < to. 

For t > to, the kinetic energy of the particles is 
too small to provide the bond breaking energy, so that 
larger and larger clusters form- We call this regime 
the aggregation regime and analyse the properties of 
the aggregates. For not too long times and sufficiently 
small volume fractions, we observe fiocculation charac- 
terized by nonoverlapping, weakly interacting clusters- 
The fractal dimension of the aggregates is approximately 
Df = 2- The cluster size distribution follows a sim- 
ple scaling form, N.m{t) ~ m~'^ f [m / fh{t)) , which has 
been applied successfully to different aggregation mod- 
els before. The increase of the typical cluster size m(t) 
can be undestood by a simple scaling analysis: Assum- 
ing that clusters irreversibly stick together when they 
hit upon each other and that their radius r grows with 
the number of particles m like r^^ ~ m, yields a cluster 
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growth m ~ f^^'/ (^-^^ 2D+2) _ r^j^jg gcaling relation shows 
good agreement with the simulation for fractal dimension 
Df = 2. 

At larger times, a spanning cluster forms, and a gela- 
tion transition is observed for all finite volume fractions. 
At the gelation transition a spanning cluster coexists 
with many small ones, wheras at very long times almost 
all particles are connected to one large cluster. On the 
largest length scales the final cluster is no longer a fractal 
but compact, as one would expect for a spanning cluster 
in the percolating phase. On smaller length scales, how- 
ever, we find fractal structures with JDj « 2. The range 
where a nontrivial fractal dimension can be observed in- 
creases with decreasing density as | logfj)]- 

Even on the longest time scales, the temperature con- 
tinues to decay. In this regime the limiting process is 
the breaking of a bond. The probability for this process 
becomes exponentially small Pbb ^ \/ /S.E/T ^-^^1^^ as 
the temperature goes to zero. Hence the cooling law for 



high temperatures is replaced by T '--^ e '^^/^ in very 
good agreement with the data. 

Several extensions of our work might be interesting. So 
far we have completely neglected all inelasticities except 
for the bond rupture. One expects the collisions at the 
hard core to be dissipative as they are in dry granular 
media. In the simplest model these could be described 
by normal restitution. Furthermore real wet grains expe- 
rience frictional forces, coupling translational and rota- 
tional motion of the grains|23|. We are not aware of any 
such studies for wet granulates. 
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